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ABSTRACT 

Timing pulses of pulsars has proved to be a most powerful technique useful to a 
host of research areas in astronomy and physics. Importantly, the precision of this 
timing is not only affected by radiometer noise, but also by intrinsic pulse shape 
changes, interstellar medium (ISM) evolution, instrumental distortions, etc. In this 
paper we review the known causes of pulse shape variations and assess their effect on 
the precision and accuracy of a single measurement of pulse arrival time with current 
instrumentation. Throughout this analysis we use the brightest and most precisely 
timed millisecond pulsar (MSP), PSR J0437— 4715, as a case study, and develop a set 
of diagnostic tools to evaluate profile stability in timing observations. We conclude 
that most causes of distortion can be either corrected by state-of-the-art techniques 
or taken into account in the estimation of time-of- arrival (TOA) uncertainties. The 
advent of a new generation of radio telescopes (e.g. the Square Kilometre Array, SKA), 
and their increase in collecting area has sparked speculation about the timing precision 
achievable through increases in gain. Based on our analysis of current data, we predict 
that for normal-brightness MSPs a TOA precision of between 80 and 230 ns can be 
achieved at 1.4 GHz with 10- minute integrations by the SKA. The actual rms timing 
residuals for each pulsar will be approximately at the same level only if all the other 
influences on timing precision (e.g. ISM, spin noise) are either corrected, modelled, or 
negligible. 
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1 INTRODUCTION 

Pulsars are stable and rapidly rotating radio sources. This 
stability ( which rivals the stability of atomic clocks on 
Earth, see lHobbs et alJl201CT ) can be exploited through pul- 
sar timing. In brief, pulsar timing works as follows: after the 
telescope surface focusses the radio signal at the receiver 
feed, the signal is amplified, sampled and digitised. Sub- 
sequently the frequency-dependent dispersion delay caused 
by the ionised interstellar medium (ISM) is removed in a 
process called "de-dispersion". Then a high signal-to-noise 
ratio (S/N) pulse profile is obtained by averaging hundreds 
or thousands of subsequent pulses in a step named "fold- 
ing". Finally this high S/N average profile is matched with 
an independently obtained standard profile (or an analytic 
template based on the data itself), in order to derive the 
time-of-arrival (TOA) of the integrated profile. 

Precise monitoring of these TOAs allows the investi- 
gation of many aspects of the pulsar, the binary system 



it inhabits and anything affecting the radio wave propaga- 
tion. For example, studies of strong gravitational field ef- 
fects in some binary pulsar systems have in the past en- 
able d some of the most constraining tests of genera l relativ- 
ity JTavlor fc Weisberd [l989l; 1 Kramer et aljfeoodY Recent 
prog ress in software (e.g . iHobbs et al.l 120091 ) and hardware 
(e.g. IVerbiest et alj|20ich has led to predictions that within 
a decade, timing of a group of millisecond pulsars (MSPs) 
will allow various characteristics of a backgrou nd of gravi- 
tatio n al waves to be determined directly (e.g. Ijenet et al.l 
12009 ; iLee et al.l |2008| ). Furthermore, the next generation 
of radio telescopes (such as the Square Kilometre Array, 
SKA, and the Five hundred metre Aperture Spherical Tele- 
scope, FAST) will have the best chance of detecting the firs t 
pulsar-black- hole binary system (e.g. ICordes et al.l I2004T )- 
Timing of such a system could determine the spin and 
quadrupole moment of a black hole and then allow direct 
tests of the Cosmic Censorship Conjecture and the no-hair 
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theorem IjKramer et al.l 12004). Pulsar timing at very high 
precision is required to achieve the aforementioned scientific 
goals, and it is more readily achieved with MSPs because 
of their short spin periods and highly stable average pulse 
shapes. Currently, several MSPs have already been timed 
at precisions down to a few hundred nanosecon ds over time 
spans of a decade or more (|Verbiest et al.ll2009T ). 

Although the integrated profiles of MSPs appear stable 
over time scales of years, there are a variety of effects that 
can affect the shape of an integrated profile on short time 
scales: multi-path propagation in the turbulent ISM, pulse 
jitter, data processing artefacts and improper calibration, 
for example. Profile variations from these effects may only 
change the pulse shape at low levels, but will cause the sub- 
sequent TOA calculation to be less accurate and precise than 
what is expected if only radiometer noise were contributing 
to the uncertainty. This will complicate timing with the next 
generation of radio telescopes since in these cases the tim- 
ing will be limited by factors other than merely telescope 
sensitivity. 

In order to investigate the level at which short-term 
instabilities in pulse shape may affect pulsar timing with 
this new generation of telescopes, we present an analy- 
sis on PSR J0 4 37— 4 715. This pulsar was discovered by 
Ijohnston et al.l (1 19931 ) and is the nearest and brightest 
MSP known, resulting in outstanding timing precision 
that has already led to a variety of inte resting results 
(|van Straten et al.|[200H ; IVerbiest et al.ll2008l ). Furthermore, 
the TOA precision o f PSR J0437-4715 o btained by current 
instruments (see e.g. IVerbiest et alj|2010] ) is already compa- 
rable to the precision future telescopes may expect to obtain 
on other less bright MSPs (see Section [5]), making it a per- 
fect target for investigations of the pulsar timing potential 
of future telescopes. 

The structure of this paper is as follows. First we de- 
scribe the observations and data preprocessing in Section [2] 
Some statistical tools are introduced in Section [3] Next we 
review the possible effects involved in profile distortion and 
present the results of data reduction in Section [4] We con- 
clude with an overview of our main findings, prospects for 
the precision timing with the next generation of radio tele- 
scopes, and a brief discussion of future research in Section [5] 



Table 1. Basic features of the selected datasets. Note that the 
2005-07-24 dataset was used only to create a timing template and 
a receiver model. 



Date 


Receiver 


Time span 


Number 


File 






(hours) 


of files 


length (s) 


2005-07-24 


MB 


8.7 


1596 


16.8 


2005-09-07 


MB 


9.0 


1500 


16.8 


2006-12-31 


H-OH 


7.4 


212 


67.1 


2007-05-06 


H-OH 


8.9 


152 


67.1 


2008-02-24 


MB 


1.0 


180 


67.1 



the data were folded in near-real time to 16.8 s for the early 
data and to 67.1 s for the later data (see Table[T|. Off-source 
observations of a pulsed noise probe at 45° to the linear feed 
probes but with otherwise identical set-up, were taken at 
regular intervals to allow for polarimetric calibration. 

For t he data processing we used the PSRchive software 
package (|Hotan et al.l 12004 ) . We removed 12.5% of each 
edge of the bandpass to avoid possible effects of aliasing and 
spectral leakage. Two models named "single axis" and "full 
reception", respectively, were used for calibration purposes 
and details will be presented in Section 14.41 Unless other- 
wise specificied in the text, we combined the polarisations 
into total power (Stokes I) and the power across the remain- 
ing 96 frequency channels. Through the following analysis, 
TOAs and their uncertainties were d etermined th rough the 
standard cross-correlation approach (|Tavlorlll992l ). with the 
fully integrated 2005-07-24 profiles (one for each observing 
band), unless otherwise s tated. Where needed , we used the 
timing model derived bv lVerbiest et al.l i|2008l ) without fit- 
ting for any parameters. 



3 STATISTICAL TOOLS 

In order to evaluate any effects on profile shape, we first 
introduce the concepts of effective pulse number and pulse 
sharpness below and then briefly illustrate their behaviour 
with pulse S/N and TOA measurement uncertainty. 



2 OBSERVATIONS 

The data used in this paper consist of five long observa- 
tions of PSR J0437-4715, taken between June 2005 and 
March 2008 at the Parkes radio telescope. Observations 
were take n with the C altech-Parkes-Swinburne Recorder 2 
(CPSR2; iHotanet al]|2006l ). The CPSR2 is a 2-bit base- 
band recorder that performs on-line coherent dedispersion 
and records two 64-MHz wide observing bands simultane- 
ously. For the data used in this paper, these bands are cen- 
tred at observing frequencies of 1341 and 1405 MHz. It also 
effectively removes RFI online by monitoring the total power 
on /is timescales and does not record any data whenever the 
power levels deviate significantly from a Gaussian. On two 
of the five days the data were taken with the H-OH receiver, 
on the remaining three da ys the central beam of the 20 cm 
multibeam (MB) receiver (|Stavelev-Smith et al.lll996l ) was 
used, as listed in Table [T] During each day of observations 



3.1 Effective pulse number 

As pulsars are weak radio sources and individual pulses 
are often not detectable, the signal needs to be folded at 
the rotation period in order to obtain profiles with suffi- 
ciently high S/N to derive precise TOAs. It is useful to check 
whether this procedure is as effective as expected. Theoreti- 
cally, the signal is expected to increase linearly with integra- 
tion length, while the root-mean-square (RMS) of the noise 
increases according to a square-root law. Consequently, the 
corresponding improvement in S/N is expected to be pro- 
portional to the square-root of the number of pulses. Given 
TV profiles with peak amplitudes of Ai and noise RMSs of en 
(i — 1, . . . iV), the single pulse S/N is: 



(S/N) ; = 



A, 



(1) 



and the S/N of a folded profile is: 
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S/N = 



I> 



E*? 



(2) 



If the profiles are identical (Ai — . . . = An; <ti = . . . = <tn), 
we have S/N oc viV. Practically, however, effects like in- 
trinsic flux variations, scintillation and system temperature 
variations cause the S/Ns of profiles with identical integra- 
tion times to differ. This causes deviations from the scaling 
rule in the processing of real data. Therefore, we define the 
effective number of pulses as: 



7V 8i 



( S/N 



VS/N n 



(3) 



where n is the number of pulses within an individual inte- 
gration, S/N is calculated from Eq. p). and S/N mcan is the 
averaged S/N for all integrations. Effectively, iV e fc is a nor- 
malised pulse number, which corrects for the varying S/N of 
individual pulsar pulses. Consequently, the measured S/N 
of averaged profiles should scale linearly with the calculated 
\/N c f c , regardless of the brightness variations of the pulses 
involved. 



3.2 Sharpness and calculation of TO A precision 

The main goals of template matching are to calculate the 
equivalent TOA of the average profile with respect to a fidu- 
cial phase provided by the template profile; and to evaluate 
the corresponding uncertainty caused by the additive noise 
of the profile. It can be carried out both in the time-domain 
and (more commonly) in the frequency-domain, by cross- 
correlating the target profile with a high S/N standard that 
was either o btained at a different observing time or created 
analytically (|Tavlorlll9920 . Theoretically, the uncertainty of 
TOA measurements, induced by the white noise of th e pro- 
file should be in the form of (|Downs fc Reichlevl 1 19831 1: 



fi ■ S/Ni 
Here S/Ni is the equivalent single pulse S/N and 



fi: 



[U'(t)] 2 dt 



(4) 



(•») 



is the pulse sharpness parameter, with U(i) the peak- 
normalised pulse waveform and 



{n(t)n(t + r)) 



(It 



(6) 



is the noise de-correlation time scale where n{i) is the noise 
function. It can be seen that the sharpness parameter ba- 
sically relates the intrinsic profile shape to the precision of 
the TOA. 

The template matching technique produces the same 
result as predicted in Eq. ([4]) , but only if there is no profile 
shape difference between the template and the observation, 
which is an ideal assumption in practice. If profile distor- 
tion (as can be introduced by any of the effects discussed 
in Section [4]) does occur, the calculated uncertainty does 
not turn out to be as good as expected for high-S/N pro- 
files. Fig. [T] shows a simulated example of this case. Here a 



< 
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Figure 1. Simulation showing the difference between ideal tem- 
plate matching (solid line) and two cases with profile distortion 
(dashed line: profile width increased by 0.5%; dotted line: width 
decreased by 0.5%). 




Figure 2. Flow chart of the stages involved in the information 
transfer of pulsar signals. The corresponding effects that can lead 
to profile shape changes are identified for each stage. 



template is created as a Gaussian and fake observation pro- 
files are formed by adding white noise to the template after 
being broadened or narrowed by 0.5%. It is clear that the 
calculated TOA errors begin to deviate from the predicted 
uncertainty once the S/N rises to values beyond 1000. Note 
that the deviations are seen to be roughly equal for both the 
broadened and narrowed cases, which indicates that it is not 
the absolute pulse shape of the observation determining the 
reliability of its TOA uncertainty, but the relative difference 
between the observation and the noise-free template, stress- 
ing the importance of reliable template profi les. The differ- 
ent re sults in the S/N— a graph s presented bv lVerbiest et al.l 
(2010) and lHobbs et al.l (2009) already suggested this to be 
the case. As their data show, even though the integrated 
profile for PSR J0437— 4715 is supposed to be intrinsically 
stable at the folding time scale, this phenomenon still ap- 
pears at high S/N. 



4 ISSUES AFFECTING PROFILE STABILITY 

In Fig. [2j we summarise the propagation path of pulsar tim- 
ing data and identify for each stage the phenomena that 
can affect TOAs and their precision. Some of these are in- 
strumental and correctable while others induce a natural 
limit to timing precision. Detailed discussions of our current 
knowledge about these issues together with more in-depth 
investigations based on our data are presented below. 
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4.1 Dispersion changes 

4.1.1 Theory 

The dispersion measure (DM) is defined as the integrated 
free electron density along the line of sight and introduces 
a delay betw een two observing frequen cies v\ and ^2, in the 
form of (e.g. iLorimer fc Kramerl l2005) 

At ~ 4.15 x 10 6 MHz 2 pc _1 cm 3 ms x (u{ 2 - v^ 2 ) x DM. (7) 

There are two ways to correct for this dispersive delay in pul- 
sar timing data. One way is to use a filterbank to split the 
observing bandwidth up into a finite number of frequency 
channels, which are subsequently dedispersed with respect 
to each other. This is called "incoherent dedispersion" be- 
cause the dispersion is performed post detection when the 
phase information is no longer available. The alternative ap- 
proach of "coherent dedispersion" performs a deconvolution 
on the Fourier trans form of the data stream, w ithout loss of 
frequency resolution IjHankins fc RickettHl975T ). The coher- 
ently dedispersed data is subsequently stored with limited 
frequency resolution, so any error in the DM value used in 
the dedispersion, will corrupt the pulse profile permanently. 
Because of turbulent motion in the ISM and the relative 
motions of the pulsar, the Earth and the ISM, the integrated 
electron density between the pulsar and Earth is continu- 
ously changing, implying that ideally the DM value used in 
on-line coherent dedispersion would be regularly updated to 
remain close to the current value. Regular updates of the DM 
values are possible with (near) sim ultaneous multi- frequency 
or wide-bandwidth observations (|You et al.l 120071 ). but the 
accuracy of determination by this method is limited by the 
system sensitivity and complicated by frequency dependent 
evolution of the profile shape. 



4-1.2 Discussion 

For the PSR J0437— 4715 data used for this paper, given 
1024 bins across the profile and a DM of 2.644 cm _3 pc, the 
smearing time (see Eq. [7J) within a 0.5 MHz wide channel 
is approximately 0.75 bins. Note that in reality the DM 
variations of PSR J0437— 471 5 on long timescales are below 
10~ 3 cm _3 pc (|You et al.l2007T ). any shape distortion induced 
by this amount of DM deviation in coherent dedispersion 
would not be detectable. Therefore, unless the DM varia- 
tion becomes significantly larger than previously observed, 
this effect does not affect our current TOA precision. Also 
note that, since this type of distortion increases the TOA 
uncertainty by broadening the profile, the uncertainty will 
still scale following the radiometer equation for future tele- 
scopes. 



4.2 ISM Scattering and scintillation 

4.2.1 Theory 

As radio pulses travel through the ionised ISM, multi-path 
scattering can cause both constructive and destructive inter- 
ference, observed in the detected signal as apparent bright- 
ening and dimming of the pulsar signal. This effect is de- 
pendent both on observing frequency and on time because 
of the relative motion of the pulsar, the Earth and the tur- 
bulent ISM. When summing frequency channels of an obser- 



vation that is affected by such scintillation, not all frequency 
channels will contribute equally to the final profile, but ef- 
fectively a brightness-dependent weighting scheme will be 
used. In the case where the pulse profile shape varies con- 
siderably across the observing bandwidth, such weighting 
will change the resulting pulse shape as scintles move across 
the observed bandwidth. 

Additionally, multi-path scattering will cause different 
delays for signals with different path lengths and thereby ef- 
fectively broaden the observed profile. For low DM sources 
observed at high frequency, the mean of this profile broaden- 
ing is non-zero but usually not significan t compared with the 
limit ed time resolution of the backend IjCordes fc Shannon! 
[2010j). Still, the change of the pulse broadening function 
(PBF) associated with either fast, stochastic variations in 
the diffractive delay, or long-term evolution of the refrac- 
tion angle, will result in instability of the profile. The pulse 
broadening will be significant ly larger f or high DM objects 
observed at low frequency (|Cordes fc Shannon! 120101 ). The 
broadening timescale can be investigated by assuming dif- 
ferent scattering models, based on which it is also possi- 
ble to reveal the intrinsic profile shape by de-convolving 
the pulse waveform with a theoretical broadening func tion 
l|Williamsonlll973l ; iBhat et al.ll2004l ; IWalker et aLll2008l ). 



4-2.2 Data investigation 

The low DM of PSR J0437-4715 (2.644 cm _3 pc) corre- 
sponds to a short broadening timescale r s (< 10 -5 ms) and 
a large decorrelation fre quency scale (estimated by A/ ~ 
1/(2t tt s ) « 0.4 GHz, see ICordes fc Lazicll200ll ; IBhat et al.l 
2004). Consequently, provided the effective 48 MHz band- 
widtljj of our data, the shape change on short timescales by 
the PBF variation should be negligible here. 

The influence on pulse shape, by the combination of 
scintillation and profile shape evolution over frequency, can 
be investigated by dividing the entire bandwidth into sub- 
bands and carrying out individual template matching. The 
templates for the sub-bands are produced with the same 
bandwidth and central frequency as the averaged profiles. 
Fig. [3] presents an example from the 2005-09-07 dataset. 
Here the effective bandwidth of 48 MHz is divided into 
six bands of 8 MHz each. The S/N is calculated from off- 
pulse bins which show Gaussian statistics. The plot shows 
that all of the matching processes have identically behaved 
S/N— cttoa curves and the same final otoa of about 22 ns. 
Statistically, the six TOAs with an uncertainty of 22 ns each, 
are equivalent to a single one of uncertainty 22/V6 ~ 9ns. 
This is the TOA precision we obtain from template match- 
ing the fully frequency-integrated observation. It means that 
across the 48 MHz bandwidth of this dataset, the TOA pre- 
cision does not benefit from conducting sub-band template 
matching. The effect could become significant when the ob- 
serving bandwidth is comparable to or larger than the scin- 
tillation frequency scale. Still, as the effect can be dealt with 
through the application of frequency-dependent template 



As described in Section|2] the total bandwidth is 64 MHz but on 
either side of the bandpass 12.5% was removed, leaving 48 MHz 
of effective bandwidth. 
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Figure 3. Sub-band template matching for the 2005-09-07 
dataset at 1405 MHz. Curves correspond to different bands and 
the double dashed line is for the entire bandwidth combined. The 
sub-band curves yield a set of TOAs which are equivalent to a 
single TOA obtained from the full-band curve. 



dspsfQ (|van Straten fc Bailesl l201ll ). The scattered power 
was mitigated during off-line processing through applica- 
tion of the correct ion algorithm implemented in psrchive 
i|Hotan et al.l 12004 ). Given the uncorrected, mean digitized 
power a 2 in each pulse phase bin, this algorithm inverts 
Eq. (A5) of JA98 to estimate the mean undigitised power 
a 2 and the mean scattered power A via Eqs. (45) and (43) of 
JA98, respectively. The effect of correction is demonstrated 
in Fig. 13] which shows the pulse profile formed from the 
2005-07-24 dataset with and without application of the al- 
gorithm, as well as the difference between the two. The de- 
creased pulse width of the corrected profile allows higher 
timing precision. Note that the distortion would not change 
significantly once the back-end settings are stable. This 
means that the TOA precision would still scale with effective 
collecting area as described by the radiometer equation. We 
therefore conclude that this effect does not limit the current 
TOA precision, and will not limit it for future telescopes 
either, which are likely to employ digitisers with a higher 
number of digitisation levels. 



matching, we conclude that it will not be a limiting factor 
to TOA precision with either current or future telescopes. 



4.3 Signal digitisation effects 

4-3.1 Phenomena 

Shape distortion induced by instrumental effects obscures 
the true pulse shape and can therefore be expected to de- 
crease the precision of timing. Two main digitisation effects 
are pertinent to a low-bit observing system. 

The first effect is caused by the underestimation of the 
undigitised power in a system with low d ynamic range (e.g. 
only 2 bits per sample). As discussed in ljenet fc Anderson! 
i|1998T ) (hereafter JA98) the earlier arrival of pulsar emis- 
sion at the high end of the observing bandwidth causes an 
increase in undigitised power, and therefore a decrease of 
the digitised-to-undigitised power ratio at all frequencies if 
the output power level is kept constant. As a result, the 
off-pulse power will be decreased in the rest of the band. 
This effect can be avoided through dynamically setting the 
output power levels, which provides the required dynamic 
output range and therefore does not result in negative off- 
pulse dips on either side of the pulse profile. 

The second artefact is caused by quantisation errors as 
a second-order distortion, and manifests itself as an increase 
in white noise uniformly redistributed across all frequency 
channels, which is induced by the increase in pulsed power 
in one part of the band (JA98). This scattered power broad- 
ens the profile and causes additional pulse shape variations 
as a function of observing frequency, which decreases the 
achievable TOA precision. 



4-3.2 Correction 

In this paper, all CPSR2 data presented were corrected for 
the low dynamic range artefact during on-line processing, by 
the dynamic output level setting algorithm implemented in 



4.4 Polarimetric calibration imperfections 

4.4.I Theory 

When a fixed-linear-feed, alt-azimuth radio telescope tracks 
a polarised source across the sky, the feed will rotate with 
respect to the plane of polarisation by the parallactic angle 
(q), which is defined as the angle between the object-zenith 
great circle and the hour circle. The change of q combined 
with the instrumental response, will result in a variation of 
the observed Stokes parameters with time. Polarisation cal- 
ibration, the aim of which is to reveal the intrinsic profile, 
will correct this time-dependent variation, but the correction 
will only be partial if any non-orthogonality of the recep- 
tors is not fully modelled. In this case, a difference between 
profiles at different q will be seen even after calibration. In 
practice, a "single axis" model considering only the differen- 
tials in gain and phase for the two linear pola risation probes 
is usually applied (e.g. lStinebring et al.lll984h . The most re- 
ce nt (here men t ioned as "full reception") model, described 
by Ivan Stratenl (J2004T) . solves the matrix description of the 
polarisation measurement equations, accounting for differ- 
ential gains and phases, as well as for coupling and leakage 
effects between the receiver feeds. 



4-4-2 Calibration 

For comparison, here we used both the "single axis" and 
the "full reception" model (constructed from the 2005-07- 
24 dataset) for calibration of the 2005-09-07 dataset. Fig. [5] 
(a)-(b) shows two hour-long integrated profiles formed from 
the 2005-09-07 dataset, covering a different range of par- 
allactic angles and calibrated according to the single axis 
model. The large differences in linear polarisation and posi- 
tion angle demonstrate the imperfection of this calibration 
model. The difference between the total intensity profile at 
the two different values of q is shown in subplot (d) of Fig. [5] 
The same profile as in subplot (a), but calibrated with the 



http://dspsr.sourceforge.net/ 
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Figure 4. Pulse profiles at 1405 MHz formed from the 2005-07-24 dataset, demonstrating the effect of digitisation scattered power on 
pulse shape. The left-hand plot shows the zoomed-in pulse profile before (solid line) and after (dashed line) application of the correction 
algorithm, while the right-hand plot presents the corresponding difference in this phase range, clearly showing the broader shape of the 
uncorrected profile. 



full reception model, is shown in subplot (c) - this profile 
is identical for both observing times, as the difference plot 
(e) shows. Furthermore, the profiles calibrated with the full 
reception mode l compare well with th e previously published 
polarimetry of iNavarro et al.l (| 19971 ). regardless of q. The 
remaining difference is clearly far less than 2% of the total 
intensity, and is therefore below the uncertainty l evel of the 
calibr ation, as quantified through simulations bv lOrd et al.l 
(2004). Further simulation shows that the profile difference 
in plot (e) will induce TOA errors of less than 30 ns. It is 
therefore clear that the full reception model removes all po- 
larisation calibration artefacts down to the level of our cur- 
rent TOA precision on PSR J0437-4715. We hence conclude 
that polarimetric calibration does not limit the current TOA 
accuracy above 30 ns. Though the application of these cal- 
ibration schemes to future interferometers (such as SKA) 
remains to be solved, such accuracy will be achievable for 
future single-dish telescopes (such as FAST). 



4.5 Interference and unknown observing system 
instabilities 

4-5.1 Diagnostic theory 

Besides the expected effects described above, there are some 
unpredictable effects that may also affect the data quality. 
Radio frequency interference (RFI) and instrumental failure 
are the most important two. Specifically in the case of an 
observing system with only four digitisation levels (i.e. a 2- 
bit system), any excess in power (as potentially caused by 
RFI) or a temporary non-linear response in the system, can 
be expected to affect the pulse shape. 

In order to keep track of any such occurrences, the 
statistics of the digitised data can be compared to those 
expected from theory, in the following way. First, the digi- 
tised data are divided into consecutive segments of L sam- 
ples and, for each segment, the number of low-voltage states 
M is counted. The digital signal processing software that is 
used to process the 2-bit data, maintains a histogram of oc- 
currences of M that is ar chived with the pulsar data for later 
use as a diagnostic tool l|van Straten fc Bailesll201ll ). When 



the voltage input to the digitiser is normally distributed, the 
ratio <£ = M/L has a bin omial distribution as in Eq. (A6) of 
iJenet fc Andersonl (|l998t ). The difference between this the- 
oretical expectation and the recorded histogram of M pro- 
vides a measure of 1) the degree to which the input signal 
deviates from a normal distribution, and 2) the degree to 
which the sampling thresholds diverge from optimality. This 
difference, called the 2-bit distortion, is given by 



D =12 [P(M/L)-U(M)] 2 



(8) 



where V($) is the expected binomial distribution and 'H(A'I) 
is the recorded distribution of M. Separate histograms of 
M are maintained for each polarisation, and the reported 
distortion is simply the sum of the distortion in each polar- 
isation. 



4-5.2 Data analysis 

The sharpness parameter can be used to examine profile 
stability over short lengths of observing time. In Fig. [6] 
the variations of sharpness and 2-bit distortion (defined in 
Eq. El) during the observing runs of 2005-09-07 and 2006- 
12-31 are shown side by side. Note that the S/N of the 
individual integrations is typically not sufficient to apply 
Eq. (J5J directly. The sharpness is therefore indirectly ob- 
tained through Eq. Q. In doing so, the white noise un- 
certainty is estimated from cross-correlation with the single 
integration and the S/N is measured by calculation of the 
off-peak RMS. As the profile has a large on-pulse duty cycle, 
only ~ 15% of the pulse period can be used to calculate the 
noise RMS, which results in a ~ 10% uncertainty of both 
the measured S/N and the sharpness estimates. The 2-bit 
distortion is calculated by following the stages laid out in 
Section gXT] 

As shown in Fig.[6]the 2- bit distortion levels were consis- 
tently low during the 2005-09-07 observation and the sharp- 
ness was close to normally distributed. The observation of 
2006-12-31, however, shows significantly larger levels of 2-bit 
distortion as well as correlated, non-Gaussian variations in 
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Figure 5. Comparison of polarimetric calibration models. Shown are two hour-long observations from the 2005-09-07 dataset at 1405 MHz 
at different parallactic angles and calibrated with the single axis (subplots a and b) and full reception models (subplot c, identical for 
both hours of observations). Here the solid, dashed, and dotted line correspond to the total intensity, linear polarisation, and circular 
polarisation profiles, individually. The difference of the total intensity profiles for the two hours of observations are shown in subplot (d) 
for the single-axis calibration scheme and in subplot (e) for the full reception model calibration. 



the sharpness parameter. Furthermore, the 2-bit distortion 
has a high degree of correlation with the sharpness levels. 
The correlation coefficient of the two time series is —0.73, 
suggesting that the profile distortion is related to the digi- 
tisation. 

To illustrate the effect these digitisation-induced shape 
changes may have on TOA precision, we provide the S/N — 
ctoa plot for the data analysed in Fig. [7] Clearly the 2005- 
09-07 data behave as expected: they follow the theoretical 
inverse relationship and worsen slightly for S/N> 1000. This 
worsening is caused by noise in the template profile, which 
was constructed from the 2005-07-24 dataset. To demon- 
strate this, we created the S/N — cttoa curve for simulated 
data, based on a (simulated) template profile with a S/N 
identical to that of the 2005-07-24 standard profile. This sim- 
ulated result is shown as the dotted line in Fig.[7]and follows 
the 2005-09-07 curve well. Ideally, therefore, a noise-free an- 
alytic template profile would be used (as in iKramer et al.l 
Il999l ). but the small-scale features present in the profile of 
PSR J0437— 4715 require advanced modelling that goes be- 



(otoa) is several factors higher than for the 2005-09-07 data 
and limits the calculated precision to the ~ 50 ns level. 



fiarntne 
Unci). 



flattening of the otoa-S/N curve in IVerbiest et al 

The S/N — otoa curve for the 2006-12-31 data displays 
much larger deviations: for equal S/N its TOA uncertainty 



The other datasets in our investigation (2005-07- 
24, 2007-05-06 and 2008-02-24) also yield well-behaved 
S/N— otoa curves like that of 2005-09-07, which shows that 
the majority of our data is not affected by this digitisation 
artefact. It must be noted that the only objective statistic 
that we found to clearly identify the 2006-12-31 data as cor- 
rupted, was the 2-bit distortion statistic defined by Eq. 1(5} 
and that this corruption cannot be corrected post-detection. 
This demonstrates the need to investigate the digitisation 
distortion in any 2-bit pulsar timing data and to exclude 
any distorted data from future pulsar timing analyses. We 
conclude that these statistics should be used to identify and 
discard corrupted data. We also note that for systems with 
more digitisation levels, this type of distortion will not affect 
TOA precision as much (or at all). It is therefore likely that 
timing with more state-of-the-art systems on both present 
and future telescopes, will not be limited by these effects. 
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Figure 6. Sharpness B (top) and 2-bit distortion D (bottom) for the observations of 2005-09-07 (left-hand plots) and 2006-12-31 (right- 
hand plots). Results of the observing band centred at 1405 MHz are shown with solid lines; those of the 1341 MHz observing band with 
dashed lines. The shape variations on 2006-12-31 identified by the changes in sharpness are clearly correlated with changes in 2-bit 
distortion, suggesting these variations are instrumental rather than intrinsic to the pulsar. The sharpness variations of the 2005-09-07 
data are close to Gaussian as their histogram (inset) shows. 
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4.6 Pulse jitter 

4-6.1 Theoretical model 

The phases of single pulses vary around an expected, av- 
erage ph ase. For some pulsars these phase variations seem 
random (|Ekers fc M offct 1968) and are called "jitter", while 
for others there are clearly s ystematic drifts called "drift- 
ing sub-pulses" (jSutton et al.Nl9700 . These phase variations 
slightly broaden the pulse shape and induce additional ar- 
rival phase fluctuation of integrated profiles. For the case 
of pulse jitter, assuming no modulation of the single pulse 
shape and intensity, the jitter-induc ed TOA scatter can be 
written as (|Cordes fc Shannonll2010h : 



Figure 7. S/N-ct t oa relations for both real data and simula- 
tions. The solid line represents the result from the 2005-09-07 
dataset, the long-dashed line from the 2006-12-31 dataset. The 
dotted line shows the theoretical prediction from Eq. J4}, and 
the short-dashed line shows the relationship for simulated pro- 
files cross-correlated with a noisy standard. All these curves are 
for the observing band centred at 1405 MHz. 



f J = 



f! j u(t)t 2 dt 



,1/2 



NjU(t)dt 



(9) 



where N is the number of integrated pulses (assuming no 
systematic drifting), /j is the width of the probability den- 
sity function (PDF) of the phase jitter in units of pulse 
width, and U(t) is the normalised pulse waveform. Note 
that this uncertainty is attributed to emission stability of 
the source and is not related to the observing hardware. 
One can see that the higher the system sensitivity is, the 
more important pulse phase jitter may become. 
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Figure 8. Timing residuals of the 2007-05-06 dataset from the 
1405 MHz side-band data. The TOAs derived from 1-min inte- 
grations are clearly scattered on a scale larger than the estimated 
errors with a reduced chi-square of 8.1. The right column shows 
the histogram of the residuals. 




1000 



Figure 9. Simulated iV e f c -S/N curves for both the jitter free and 
jittered profile cases. The solid, long dashed, short dashed, and 
dotted lines correspond to jitter factors /j (see Eq. 2) of 0, 0.1, 
0.3 and 0.5, respectively. 



4-6.2 Data analysis 

Pulse jitter can be studied by a few methods and in the 
following we describe each approach in detail. 



4.6.2.1 Timing rms The first way to evaluate the im- 
pact of pulse phase jitter on timing, is to study the random 
variations of the TOAs after a timing model has been sub- 
tracted. When we consider the timing residuals for the 2007- 
05-06 dataset, we notice that the TOAs are widely scattered 
and the reduced chi-square is well above unity (shown in 
Fig. [8)l . This is so because the TOA uncertainties have been 
determined based purely on the amount of radiometer noise 
in the pulse profile, while other possible contributions of un- 
certainty such as phase jitter, timing model imperfections, 
short-term interstellar instabilities and instrumental effects 
remain unquantified. Note that a previous single pulse stud y 
has shown no evidence for pulse drifting IjJenet et al.l l 1998). 
The dataset also has passed through the 2-bit distortion test 
(mentioned in Section 2321 an d the residuals satisfy a dis- 
tribution close to Gaussian, which suggests the insignificance 
of effects such as faulty ephemeris and improper calibration 
that can induce non-white noise in timing. The deviation 
from Gaussian distribution is attributed to both small num- 
ber statistics, and to differences in measurement precision 
of the residuals caused by S/N variations in individual in- 
tegrations. Consequently, we can set up an upper limit on 
the phase jitter contribution, by assuming that jitter noise is 
dominant in the additional residual scatter (the amount not 
quantified by the TOA uncertainties from radiometre noise) 
and using Eq. ([5| to derive the parameter /j. We find that 
in this worst-case scenario, /j ~ 0.08 and <TR a d/o"j — 0.3, 
which means that the timing residuals are dominated by 
pulse phase jitter and can therefore hardly be improved by 
increasing the telescope gain, but only by extending the in- 
tegration timqj. 



3 The result is consistent for the other datasets and a detailed 
study will be provided in a subsequent paper 



4.6.2.2 A^efc-S/N relation A second approach to inves- 
tigate pulse jitter is through the N e { c — S/N relation intro- 
duced in Section T3. II A deviation from its theoretical scaling 
can be induced by inaccurate folding, non-Gaussian noise in 
the off-pulse baseline or, more interestingly, by pulse jitter. 
Based on the model of Eq. ([9]), the effect of pulse phase jitter 
on the x/A^ofc-S/N scaling law has been evaluated through 
Monte-Carlo simulations, as shown in Fig. [9] Here the tem- 
plate of PSR J0437— 4715 is used as a single pulse shape, and 
for each pulse we apply a Gaussianly distributed shift based 
on a given /j value. Clearly, the inclusion of pulse phase 
jitter introduces an initial deviation that is strongly depen- 
dent on the value of /.j. After integration of a sufficiently 
large number of pulses (of the order of a few to several tens 
of pulses), the scaling law is recovered while the curves re- 
main at a lower S/N than in the jitter-free case. Fig. \W\ 
complements Fig. [9] by plotting the same parameters over 
a range of increasing integration times. While Fig. [9] shows 
the simulated effects of pulse phase jitter at the very shortest 
integration lengths that are inaccessible to our data (which 
have a minimum length of 16.8 seconds or just under 3000 
pulses), Fig. 1 101 shows the (expectedly stable) behaviour on 
longer timescales. The S/N re f used for the calculation of 
A'cfc is 150, which is roughly the mean value for all 67.1s 
integrations. The left plot underlines how the S/N does not 
scale linearly with the real number of pulses averaged while 
the right plot does follow the theoretical linear relationship 
as expected. The phenomenon shown is visible only when 
the shortest integration time can resolve of the order of ten 
pulses. These plots do not show any evidence of the effects 
simulated in Fig. [5] This is to be expected, as the simulations 
show that the jitter effect saturates after a relatively short 
integration time, which is not accessible by the (relatively 
long) integrations of the data we have. 

4.6.2.3 Sharpness variation A third diagnostic plot 
that can be used to analyse pulse phase jitter, is the fi — AT e f c 
plot. An example of this is provided in Fig. 1111 which shows 
how fi converges as more pulses are integrated. This sim- 
ulation (which was performed in a similar way to that de- 
scribed above) is for the extreme jitter-dominated case. Real 
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Figure 10. Real and effective pulse number versus S/N, demonstrating the relations between integration time and S/N. These plots are 
based on data at 1405 MHz from 2005-07-24 (solid line), 2005-09-07 (long dashed line), 2007-05-06 (short dashed line) and 2008-02-24 
(dotted line). The theoretical scaling law is clearly reproduced in the JV e f c graph but not in the true N graph. 
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Figure 11. Simulated integration process in the existence of 
pulse jitter. The profile shape estimated by the sharpness pa- 
rameter is well constrained below the 0.5 % level after folding of 
a few hundred pulses. Note that given the rotational period of 
PSR J0437— 4715, 10 4 pulses roughly correspond to an integra- 
tion time of 1 minute. 



data would be a combination of the jitter-induced exponen- 
tial convergence simulated in Fig. QT] and of Gaussian white 
noise. However, because the intrinsic shape of single pulses 
for PSR J0437— 4715 is not easily defined, we cannot com- 
pare Fie. I 111 easily with the fi values obtained from our data 
(as shown in Fig. [6}, but we do point out that fi does not 
change significantly as we integrate our data. Furthermore, 
the variation of fi shown in Fig. [6] (top left panel) is far larger 
than the variation suggested by the simulations (Fig. Ill[) . 
This suggests either that the simulated level of pulse phase 
jitter is well below the radiometer noise present in the data, 
or that there are significant variations in the shape of single 
pulses, where our model assumes there are none. 



4-6.3 Conclusion 

In summary, we have used timing measurements to de- 
rive the upper limit of /j ^ 0.08 for pulse phase jitter in 
PSR J0437— 4715. This upper limit was subsequently used 
in simulations to evaluate the impact such jitter would have 
on the pulse shape and we found that at integration times 
beyond 100 pulses (~ 0.57 s) neither the S/N of the inte- 
grated profile, nor the pulse sharpness are affected signifi- 
cantly compared to the radiometer noise. This lack of pulse 
jitter effects on the S/N and sharpness of the simulated data, 
implies either that the derived upper limit is overly conser- 
vative and that the underestimation of TOA uncertainties 
is mostly induced by causes other than pulse phase jitter, or 
that our model for pulse phase jitter is overly simplistic and 
a more complex model is required. A more detailed investi- 
gation is required to distinguish between these two scenarios 
and to quantify the pulse phase jitter in PSR J0437— 4715 
more precisely. Such an analysi s is being prepared for a sub- 
sequent paper (|Liu et al.ll201ll ). 



5 DISCUSSION AND CONSEQUENCES FOR 
FUTURE TELESCOPES 

5.1 Summary 

In this paper, we have used the brightest and most precisely 
timed MSP, PSR J0437-4715, to illustrate the most impor- 
tant phenomena known to affect the shape of pulse profiles; 
to estimate the impact of these phenomena and to evaluate 
the efficacy of mitigation schemes. By concentrating on the 
brightest MSP currently known, we are able to cast some 
light on effects that future telescopes like FAST or the SKA 
will come across as a matter of course in any standard MSP. 
We find that pulse phase jitter may be dominant in cur- 
rent short-term timing of PSR J0437— 4715, though further 
analysis is required to enable quantification of this effect. 
The effects of faulty de-dispersion are found to be small, 
negating the need for frequent updates of the dedispersion 
DM. The effects of scintillation across the bandwidth could 
not be fully investigated given the limited bandwidth of our 
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data, but usage of frequency- dependent templates can be 
expected to resolve the aforementioned problem that scintil- 
lation might cause. We present the results of the application 
of correction algorithms for the most important digitisation 
artefacts and illustrate the importance of full calibration 
modelling, as opposed to the more traditional calibration 
for differential gain and phase only. Finally, we point out 
the importance of adequate monitoring of the digitisation 
statistics for 2-bit pulsar observing systems, which is not 
commonly practised. 

We further propose a few diagnostic plots to assess the 
data quality of any pulsar timing data: 

• Time-fi curve: Any non-Gaussian variations in profile 
sharpness suggests changes in pulse shape and imply the 
data should be carefully studied before being included in 
any timing analysis (e.g. Fig. [6}. 

• S/N-otoa relation: To assess the quality of the stan- 
dard profile used (Fig. [1} . Deviations from the theoretical 
relationship or significant differences between datasets indi- 
cate data distortion (e.g. Fig. [7]). 

• S/N-cttoa for sub-bands: As in the previous point, any 
mismatch between curves of different sub-bands indicates 
loss of data quality (e.g. Fig. [3J. 

Note that the analysis is applicable not only to future tele- 
scopes, but also to existing ones such as Parkes, Arecibo, Ef- 
fels berg, and the Large E uropean Array for Pulsars (LEAP, 
see iFerdman et al.l I2010T ). if the system sensitivity is high 
enough given the brightness of the source. 




Integration time (min) 

Figure 12. Expected average-brightness MSP TOA scatter ver- 
sus integration time for different instruments. Only radiometer 
noise and pulse jitter are accounted for as influences on the mea- 
surement accuracy. 



Table 2. Key parameters used in prediction of TOA precision for 
different instruments at 1.4 GHz. 



A eB (m 2 ) 



(K) 



Parkes 
FAST 



2.2 X 10 3 
7.1 x 10 4 



SKA core 2.5 x 10 5 



23.5 
20.0 
30.0 



5.2 Timing in the new millennium 

Since PSR J0437— 4715 is two orders of magnitude brighter 
than most MSPs, current observations of this pulsar provide 
a good demonstration of the TOA measurement situation for 
most of the MSPs by the next generation of radio telescopes. 
The 21 st century radio telescope will significantly increase 
the S/N of pulsar detections and correspondingly reduce the 
uncertainty of TOA measurenrents through vast increases in 
effective collecting area. 

Fig. [T5] shows the expected MSP TOA precision that 
can be achieved by SKA and FAST. Here for Parkes and 
FAST we only consider the uncertainty induced by radiome- 
ter noise and pulse jitter at the worst-case level derived in 
Section 14.6.21 while for the prediction of SKA we calculate 
the upper limit by applying the worst-case in jitter and set 
up the lower limit by assuming no pulse jitter at all. The in- 
strumental effects causing profile distortions and phase fluc- 
tuations are neglected based on previous analyses and the 
ISM influe nce on pulse shape is assumed to be either cor- 
rected (e.g. IWalker et alj|200g ) or not significant on a short 
timescale with provided frequency and bandwidth. Concern- 
ing the properties of MSPs, we assume a mean flux density 
of 3.0 mjy at 1.4 GHz, 50 MHz bandwidth, spin period of 
5.0 ms and 5% pulse width. In addition, a 10 K sky tem- 
perature at 1.4 GHz is applied. Consequently, it is shown 
that the TOA precision for normal MSPs can be improved 
by over one order of magnitude with the next generation of 
radio telescopes, even when assuming a worst-case level for 
pulse phase jitter. The result indicates that jitter induced 
uncertainties will be considerable in future timing of MSPs. 



The instrumental par ameters used for the s e calculati ons are 
presented in Table l2l ()Schilizzi et al.l[2007l : lNanll2006h . 

We conclude that at 1.4 GHz for 10-min integrations, 
a TOA precision of between 80 and 230 ns can be expected 
in timing of normal brightness MSPs by SKA. The wide 
range of this prediction is caused by our limited knowledge 
of the pulse jitter mechanism. The future radio telescopes 
will enable deeper investigation of profile stability of single 
pulses. Presently the comparatively high levels of radiome- 
ter noise mean that only a few bright MSPs have allowed a 
systematic study of their single pulses, wher e only a subset 
of the single pulses can be clearly detected |Cognard et al.l 
1 19961: IJenet et al.lll998l : lEdwards fc Stappersll2003l ). The sig^ 
nificant improvement in instrumental sensitivity will both 
dramatically increase the number of MSPs available for sin- 
gle pulse analysis and reduce the selection effects present in 
current studies. 



5.3 Limitations and future work 

The analysis based on PSR J0437-4715 is limited by the 
selected sample. Note that the DM and scattering timescale 
(see Section 14. 2p of this pulsar are relatively low so that 
the influence of the ISM on the TOA precision within the 
provided bandwidth is expected to be negligible. As system 
sensitivity and observational bandwidth are improved, more 
high-DM sources will be timed at high precision TOA mea- 
surements. The effects of scattering and scintillation on pro- 
file shape will then become more considerable than shown 
in this paper, and will need to be accounted for by using 
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correction algorithms JBhat et alj|2004 IWalker et alj|2008l l 
and/or frequency-dependent template profiles. 
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